%% Find and Describe Balanced Growth Path
% by Jaromir Benes
%
% The SPBC.model is a BGP model: It does not have a stationary long run.
% Instead, it has two unit roots, introduced through the productivity
% process, and the general nominal price level. To deal with BGP models,
% there is absolutely no need to stationarise them. They can be worked with
% directly in their non-stationary forms.

%% Clear Workspace
%
% Clear workspace, close all graphics figures, clear command window, and
% check the IRIS version.

clear;
close all;
clc;
irisrequired 20140315;

%% Load Solved Model Object
%
% Load the solved model object built in `read_model`. Run `read_model` at
% least once before running this m-file.

load read_model.mat m;

%% Compute Two Different Points on BGP
%
% Compute two different points on the BGP corresponding to two different
% levels of productivity, `A`. The resulting steady-state levels of other
% varibles are always in constant proportion to the level of `A` (here,
% they simply double). The steady-state growth rates remain, obviously,
% unchanged. Whenever some variables are fixed in `sstate`, the steady
% state must be solved for non-recursively, i.e. with `'blocks=' false`
% <?blocksFalse?>; this is the default setting, and can be therefore
% omitted.

oo = {'tolX=',1e-16,'tolFun=',1e-16};

m1 = m;
m1.A = 2;
m1 = sstate(m1, ...
    'growth=',true,'blocks=',false,'fixLevel=','A','display=','final', ...
    'optimSet=',oo); %?blocksFalse?
chksstate(m1);

m2 = m;
m2.A = 4;
m2 = sstate(m2, ...
    'growth=',true,'blocks=',false,'fixLevel=','A','display=','final', ...
    'optimSet=',oo); %?blocksFalse?
chksstate(m2); 

disp('Productivity level and gross rate of growth')
m1.A
m2.A

disp('Output level and gross rate of growth')
m1.Y
m2.Y

disp('Real wage level')
real(m1.W) / real(m1.P) %#ok<NOPTS>
real(m2.W) / real(m2.P) %#ok<NOPTS>

%% Solve Model Around Different Points
%
% It does not matter which point on the BGP is used to solve the model.
% They give the same solution. Illustrate this fact here by comparing the
% covariance matrices of the model variables, and a shock simulation.

m1 = solve(m1);
m2 = solve(m2);

C1 = acf(m1);
C2 = acf(m2);

index = isfinite(C1);
maxabs(C1(index),C2(index))

d1 = zerodb(m1,1:20);
d1.Er(1) = 0.01;
s1 = simulate(m1,d1,1:20,'deviation',true);
s1 = dbextend(d1,s1);

d2 = zerodb(m2,1:20);
d2.Er(1) = 0.01;
s2 = simulate(m2,d2,1:20,'deviation',true);
s2 = dbextend(d2,s2);

[s1.Y,s2.Y,s1.Y-s2.Y] %#ok<NOPTS>
maxabs(s1,s2)

%% Help on IRIS Functions Used in This File
%
% Use either `help` to display help in the command window, or `idoc`
% to display help in an HTML browser window.
%
%    help model/sstate
%    help model/solve
%    help model/acf
%    help model/subsasgn
%    help model/subsref
%    help model/zerodb
%    help model/simulate
